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Abstract 

This preliminary worl^ presents a simple extension of the explicit Euler 
method in order to improve the accuracy of the solution for relative high 
time-step. The proposed method uses an explicit Euler scheme to compute 
firstly an inaccurate solution. Then, a decomposition of the solution into 
specific convergent series allows to refine the precision of the solution that 
does not depend a priori on the chosen time- step. 



^This work is distributed under CC license http://creativecoinmons.org/licenses/ 
by-nc-sa/3.0/| 



1 Introduction 



Despite its lack of accuracy that prevent it from being used as an ODE solver, the 
Euler method is known to be very easy to implement [Ij. In particular, many efficient 
methods have been successfully implemented to solve ODE in a parallel way. For 
example, [2] introduced the possibility of parallelizing the computations requested to 
solve efficiently an Euler scheme. This method uses a predictor-corrector algorithm, 
applied to an implicit Euler scheme, and allows a faster simulation, especially if 
parallel computing is available p]. Further analysis and stability proofs can be found 
e.g. in [3] [1] [5] [6]. Another parallelization methods and comparisons are also 
discussed in e.g. [7][8]|9][Tn][n][l2][l3][Tl]. 

The proposed method is based on the explicit Euler scheme for which we aim 
to improve the convergence at each time-step. Based on the explicit Euler scheme 
properties, a specific series expansion of the solution at each time-step allows to get 
a more accurate precision of the solution. 

The paper is structured as follows. Section II presents the outhne of the method. 
Section III gives a brief overview of the acceleration algorithm used to increase the 
precision of the solution. Simulations results and some possible improvements may 
be found in Section IV. Some concluding remarks may be found in Section V. 



2 Outline of the method 

Consider a partial differential equation, eventually non-linear, such as: 

f^=A{t)y{t) + B{t)u{t) te[0,Tf] (1) 

with the initial condition y{0) = yo. u and y represent respectively the input and the 
solution of (jl]), which can be a scalar or a vector. To solve ([T]) and compute an initial 
solution over [0, Ty], we use the explicit forward Euler method for which we assume 
initially that the time-step is large and ensure the stabilitjj^ The equation ([T]) is 
rewritten in the discrete domain: 

7 - ^koVko + ^koUko [^) 

whose solution yl^+i verifies for the ko + 1 time-iteration: 

yt+i = + MMVko + BkohkoUko ^ = Sho ivt'^ko) (3) 

*This condition is a priori not necessary and may be investigated in a future work about the 
stability of the proposed method. 
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The coefficient ho is the initial time-step such as the solution is calculated at each 
instant t^^ = koho, G N. and B^g are "connected" to A(t) and B(t) and are 
described by specific relationships that depend on the smoothness of A(t) and B(t). 
We call / + Afcg/io, the dynamic matrix of the explicit Euler scheme. Symbolically, 
we introduce the operator S, as the Euler iteration in order to increase the solution 
from to corresponding to the ko + 1 step, y^^ is thus considered as an initial 
condition of the ko + 1 step. 
We denote: 

^° = {Z/O, y\tl), y\tl) , ■ ■ ■ , y'im 

the set of the calculated solutions (with the time-step ho) at each instant = koho. 
Consider now a smaller time-step hi = ho/ Si < ho, with > 1, for which we obtain 
the set of the calculated solutions: 

= {yo, y'itl), y\tl), ■ ■ ■ , y\tl) , ■ ■ ■ , y\Tf)} 
at each (same) instant = kihi = kiho/Si = koho, ki > ko. 

Therefore, by induction, we can deduce that for any time-step hi = ho/ Si < with 
Si > Si-i > 1, the set of the calculated solution^ 

= {yo, y\t% y\tl), ■ ■ ■ , ■ ■ ■ , y\Tf)] (4) 

at each (same) instant = kihi = kiho/Si = hi_iki_i = hoko, ki > ki_i. 
Similarly, we denote by the set of the "true" solutions: 

^^ = {l/0,/(t?),/(t°),---/(O, •••,/(T,)} 

at each (same) instant t1^. Note that since ho is the initial time-step and thus the 
"reference" time-step, y^ is computed at the same instants t1^. 
Therefore, at each instant t^^, ko G N, the solution is described as an (infinite) series 
composed of terms that are computed / deduced only from the dynamic matrix 
(/ + Ak^ho) powers. Basically, at each instant the solution ([3| is described as 
a symbolic series (composed of the underlined terms of each \E'*,'i G N): 

^^;.o = (y'itl), y\tl), ■ ■ ■ , y\tl) • • • ) (5) 

that converges to ^ V^i^^ko)- obtain an accurate estimation of the limit fi^*™, 
we will describe the e-algorithm, which is a convergence accelerator algorithm, whose 
purpose is to compute the limit from only a few terms of VL^q. 

^In a similar way, using the £ operator, we may have : 
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Illustrative example : Consider a linear second order system Ei, described by a 
state-space representation: 



Si := 



Xi 

i2 



-1 5 
-5 -1 



Xi 
X2 



10 



a;i(0) =0,X2(0) 



(6) 



Figure [T] presents the different solutions of Si, relating to the different time-step and 
evolving in the phase-space Xi — X2. The "true solution" is the exact solution of Si. 
Starting from a time-step ho, that gives a very inaccurate solution, the solution is 
also computed considering 6i,i = 1...5. 
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Figure 1: True solution of Si plotted in the phase-space in comparison with explicit 
Euler scheme with different time-steps hi (corresponding to the "ifi series). 

Therefore, at each instant tl^,ko G N, we define the series such as (denote 
(xi X2)'^ , the state-space vector): 

_ /0.5000 1.4000 2.4900 3.4980 4.1570\ 
°~ VI. 9000 2.4600 2.5140 2.0176 1.0668^1 ° ~ 

'0.6000 1.4911 2.4295 3.1761 3.5565\ 



' 1.8150 2.2146 2.1270 1.6077 0.8168J ^"^^ 

_ '0.6256 1.5013 2.3850 3.0589 3.3806\ , _ , /o 
^ ~ ' 1.7805 2.1333 2.0185 1.5152 0.7893 / ^ ~ °/ 
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_ /0.6427 1.5030 2.3440 2.9668 3.2522\ , _ , 
^~\1.7517 2.0705 1.9404 1.4552 0.7808^ ^^^" °/ 

_ /0.6492 1.5021 2.3254 2.9281 3.2007\ 
^~V1.7391 2.0445 1.9094 1.4329 0.7799^1 ^ ~ °/ 

and the resulting series ilk^ (denote (xi 0:2)^, the state-space vector): 



_ /0.5000 0.6000 0.6256 0.6427 0.6492\ _ /0.6630 

i^fco=i - I i_9ooo 1.8150 1.7805 1.7517 1.7391 j ^ '=0=^ ~ V 1.7075 



1.4000 1.4911 1.5013 1.5030 1.5021\ /1.4963 



^fco=2 ( 2.4600 2.2146 2.1333 2.0705 2.0445/ ^^^0=2 I 1.9817 



_ /2.4900 2.4295 2.3850 2.3440 2.3254\ _ /2.2770\ 

iifeo=3 - I 2.5140 2.1270 2.0185 1.9404 1.9094) ^ '=0=^ ~ I 1-8379 j 



3.4980 3.1761 3.0589 2.9668 2.9281\ . ^^m /2.8346 



^ko-A (2.0176 1.6077 1.5152 1.4552 1.4329] ^^^0=4 I 1.3851 



4.1570 3.5565 3.3806 3.2522 3.2007\ /3.0809 



[iQQQg 0.8168 0.7893 0.7808 0.7799 J ^^^0=5 \^o.7836 

At each instant t^^, ko e N, given some terms of the series Q^g, the use of an accel- 
eration algorithm, allows to have an accurate estimation of fl^J^- We call standard 
Euler-Shanks scheme, the application of the Shanks transform to the explicit Euler 
scheme. 



To solve accurately the explicit Euler scheme with a high (stabilizing) time-step, one 
performs a first resolution using a high time-step Hq. Then, for each instant kho, k eN 
of the initial solution, some power computations of (I -\- Ak^hi) at different time-steps 
hi give a sequence for which, the limit, and therefore, a good estimation of the true 
solution, is deduced using a convergence accelerator. 
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3 Shanks transform of numerical series 



Consider {Sn)neN (also simply Sn), a real series that converges to a limit Sum for 
n > niim- To accelerate the convergence, one defines the transformation $ : M — > M, 
such as Tn = ^{Sn), whose limit is Too, implies that Too — Sum = o{Sn — Sum)- In 
other words, we define the transformation $, such as the transformed series T„ from 
the series Sn reaches its limit Sum faster than Sn- 

Among the different methods that has been established to accelerate the conver- 
gence of a series [15] [IE], we consider the e-algorithm [T7] [18] [I9] based on the 
Shanks transform [211] [21] [22]. 



3.1 The ^-algorithm 

Consider a series Sn, for which we aim to estimate the limit Sum- From a few terms 
of Sn, the limit Sum can be extrapolated using the following algorithm[^ 



4"^ = 5„, (neN) 
where Sn is either a scalar or a vector. 

The e-algorithm can be described as a function S** = Se(k,n, Sn) that extrapolates 
Sn and gives a new series S** that converges faster to Sum (the limit Sum should be 
reached after a few terms of S**), where Sn is the series for which we aim to estimate 
the limit Sum] k and n are respectively the index of Se {k is also the order of the 
associated Shanks transform) and the (initial) index of Sn- 

Table [T] gives the number of terms of Sn and the corresponding number of computa- 
tions that are requested in order to compute Se{k,n, Sn) - Only even k number are 
considered and any index n from Sn may be considered {{k,n) G N). 



Table 1: Table of computations of Se{k,n, Sn)- 

k nb. of operations Sn Sn+l Sn+2 Sn+3 Sn+i Sn+5 Sn+G 



1 1 - 

2 6 2 3 

4 35 3 10 

6 204 4 21 

8 1189 5 36 



1 












14 


7 


1 








55 


70 


42 


11 


1 




140 


301 


363 


242 


86 


15 



From this Table, we can deduce that the estimation of the limit Sum is more 
accurate when k and n are sufficiently high (depending on the convergence rate of 

"''To compute the inverse in the vectorial case, considering a vector U, we have : U^^ — jytjtyj- 
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Sn)- However, higher index involve much more computations that may decrease 
the overall efficiency of the method. As a result, a compromise has to be found 
between k and n in order to get the most accurate Sum with a few terms of Sn- Two 
improvements could be considered: 

• From the index n, instead of considering the consecutive terms Sn, 5'„,+i, Sn+2, 
■ ■ ■ we define the map Q : Sn ^ S^ such as the application 9(5'„) modifies the 
series Sn in order to rearrange "pertinent" terms. For example, given the series: 

Sn = {'S'n, Sn+li Sn+2, Sn+3, ' ' ' } 

a linear map 0(5'„) gives: 

Sl^ = {Sn, Sn+A, Sn+8, Sn+12, ' ' '} (§) 

and a non-linear map would give for example: 



S'^ — {5'Lexp(n)J, Slexp{n+1)\, 5'Lexp(n+2)J , 5'Lexp(n+3)J ; " " " } (9) 

where [.J is the floor function defined by \_x\ = max {m G Z | m < x}. 

The purpose of the linear map is to distribute uniformly the terms of the series 
Sn, whereas the non-linear map may emphasize the high index terms of Sn- 

• To use more efficiently the S^ function and therefore to get a more accurate 
estimation of Sum, one can consider repeated actions of S^ on Sn 



3.2 Application to the explicit Euler scheme 

To provide a "generic" form of the explicit Euler induction, we assume that A, B and 
u are constants using the "initial" time-step Hq- From the ffist recursive terms of (|3|: 

y'i = {I + h^A)y^^ + hoBu 
yl = {I + hoA)yf + hoBu 
yl = {I + hoA)y^ + hoBu 
yl+i = {1 + hoA)yl + hoBu 
we obtain the following definition of the explicit Euler scheme: 

fco-l 

yl^, = (J + Aho)'° + ^ (/ + A hoY hoBu, heN (10) 

i=0 
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According to the Section|2| to describe the exphcit Euler scheme with lower time-step, 
we introduce the parameter that divides the time-step ho such as hi = ho/Si, and 
the time-iteration /cj, as a subdivision of /cq- Therefore, on the interval fcj G [ko, ko + l], 



the equation (10) is therefore rewritten: 



= (/ + A^y^ yl + + jJj:Bu, h e N (11) 

The notation y].-j^i refers to the fact that the solution is computed at the ki + 1 
th time-iteration with the time-step h^/di (the referenced time-step is ho). In other 
words, we deduce and retrieve the fact that if ki = 6i, then t^. = kiho/Si = h^. It 
follows that the corresponding fi^g (defined by (IS])) vectors are defined by: 



^;^o = ivk UMo)^ yl \sMo)^ ■■■^yk IsAy^.) ■ ■ ■ ) (12) 

Since the terms of flko are directly computed from the parameters 6i,i G N, we can 
define in the same manner, the vector A^,, that is in direct correspondence with Q^g-. 



^fco 



{0,61,62, ■■■ ,Si,---) (13) 
Then, for each time-iteration ko, we define the accelerated vector Q*^^ such as : fl^^ = 

Ss{k,n,Qko)- The estimated limit of is noted r^^g*™" and is given by the last 
term(s) of Ql^. 

4 Examples and possible improvements ... 
4.1 Linear second order system 

Consider again the linear second order system Si, described by the state-space rep- 
resentation: 



- 1 ) = ( Is -1 ) ( ::: ) + ( ? J w -.(o) =o,..(o) = 1 



Instead of applying the mapping 9 on fi^g defined by (12), a more convenient way is 



to map the vector defined by (13). Therefore, starting from a linear A^g vector. 



we can consider different map as examples. For example: 

6o = I 6i+i = 6i + l 

Tables [2| [3] and |4] give details about the execution of the e-algorithm for ko 
considering different k index and different maps G. 
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Figure 2: Application of the e-algorithm scheme to the explicit Euler integration of 
Si. The big black circles represent the estimated limits f^^*™", (^o = 1...5). 



Table 2: Execution of the e-algorithm for a G linear mapping (|8j) composed of 10 
terms of flko with 6io = 400. 



k 


min5j G A 


max 6i G A 


error with the true solution (%) 


4 


1 


178 


0.154 


6 


1 


267 


0.101 


8 


1 


355 


0.073 



Table 3: Execution of the e-algorithm for a G linear mapping (|8]) composed of 10 
terms of flko with 6io = 300. 



k 


min 6i G A^^ 


max^j G Afco 


error with the true solution (%) 


4 


1 


133 


0.204 


6 


1 


200 


0.133 


8 


1 


266 


0.096 



We deduce that the choice of the map G and the index k are very important to ensure 
a good accuracy of the estimated limit. An optimization process can be therefore 
applied to define optimal G and k according to the smoothness of the system to solve. 
Nevertheless, good accuracy can be obtained for a small amount of \&j computations 
with a minimal time-step of the range of /iq/IOO, ho = 0.1 in our case. 
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Table 4: Execution of the e-algorithm for a O exponential mapping ([9]) composed of 
10 terms of ilkp with Sip = 22026. 



k 


min 6i G Afc„ 


max 6i G A/jjj 


error with the true solution (%) 


4 


1 


148 


0.555 


6 


1 


1096 


0.222 



Applying the e-algorithm to each vector fi^p, (/cq = 1...5) results in^the estimated 
limits f^fc™, (^0 = 1---5) represented by the big black circles in Figure 

4.2 Considerations about non-linear systems 

Consider a non-linear first order system ZI2, described by the OD 

S2 := + (t^-t- 3)y(t) = 3 sin(t - 0.25) (15) 

at 

Equation ( 15 ) is of the form of ([T]) for which : 



A{t)=r-t-3 S(t) = 3sin(t-0.25) 

The use of an initial high time-step ho may introduce some errors in the integration of 
A{t) and B{t), that prevent them to be defined as constant during an Euler iteration. 
To discretize A{t) and B{t), we define an interpolation function that may linearize 
([1]) at fco- A possible form is: 

where a, /3, 7 and 6 are real parameters that are adjusted according to the smoothness 
of Ako and 5^^. 

Figure [3] presents the computation of the solution of S2 with different time-steps in 
comparison with an high-resolution explicit Euler {h = 10"''), that we assume in 
this the "true solution". The accuracy of the linearization depends on the 

interpolation function and its parameters. In particular, if /3 = 5 = 0, then the 
solution is far from the "true solution". 

As a result, we apply the same Euler-Shanks scheme to the linearized equation 



Example taken from the Matlab® documentation about " ODE solvers" . 
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Figure 3: S2 computed using a high resolution exphcit Euler {h = 10 ^) in comparison 
with Explicit scheme with different time-steps. 



4.3 Pre-Adaptive scheme 

We propose a pre-adaptive scheme, that uses, at each instant t^^, the previous esti- 
mated limit f^^Q*™ to recompute the initial solution of ([l]). 

Starting from the initial condition yo, a first Euler iteration is computed such as 
Vi = ^hoiVoyUo) for the instant t^^^i- Then, several Euler iterations Sh^ using different 
time-step 5, give solutions at the same instant t1^=i for which, as a result, the vector 
ftko=i is computed. The e-algorithm applied to ^^0=1 gives the accelerated vector 
^Iq=i and thus the estimated limit 

To compute the second Euler iteration, instead of taking as the initial condition, 
we take ^ij^^^i- Computing new Euler iterations (for different Si) produce a vector 
^fco=2 at the instant t°g=2 5 which the estimated limit ^^fc^!!^ become the initial 
condition of the next Euler iterations... 
As a result, using the £ operator, Q is rewritten: 

K = {yo,ShXyo,uo), f,,(fi^;L"l, wi), ■ ■ ■ , s,X^llr,,u,,.,) ■■■} (16) 

To illustrate the application of the pre-adaptive scheme on the system Si, we 
consider a G linear mapping ([s]) composed of 10 terms of ilko with Siq = 400. Table 
|5] presents a comparison between the "standard" accelerated Euler scheme and the 
pre-adaptive accelerated Euler scheme for ko = 2. 
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Table 5: Comparison of the "standard" scheme and the pre-adaptive scheme. 



k 


min 6i e Ak^ 


max 6i e 


error (%) - 'standard' 


error (%) - pre-adaptive 


4 


1 


178 


0.154 


0.077 


6 


1 


267 


0.101 


0.050 


8 


1 


355 


0.073 


0.036 



Since the first term of each vector is already defined by the previous estimated 
limit ^ko^Il, then the following estimations are more accurate and the e-algorithm 
would require less terms (low k index) than in the "standard" case. Figure |4] presents 
the integration of the Si system with the initial time-step ho in the case of the 
"standard" scheme and the pre-adaptive scheme. "$1^ converges a priori to the true 
solution in a few iterations. 

Moreover, one can consider also the application of the repeated e-algorithm [23] 
considering higher index k for the first terms of 




Figure 4: Pre-adaptive scheme applied to Si in comparison with the 'standard' case. 



5 Conclusion 

We presented a simple derivation of the explicit Euler scheme that uses a convergence 
accelerator algorithm in order to improve the precision of the resulting solution at each 
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instant. Simulations show encouraging results and show that the proposed method 
has the following features: 

• The computations are based essentially on the power computation of the dy- 
namic matrix (/ + A^^hi) (thus, no matrices inversion). 

• Since the power computations involve several time-steps (the initial time-step 
ho, and the different division 5, requested by the e-algorithm) , these computa- 
tions may be parallelized. 

• Multidimensional ODEs (linear or non-linear) can be a priori considered. 

Future investigations include the study of the stability of the proposed algorithm, 
the improvement of the series accelerator algorithm [19] |23] and the adaptive scheme 
(e.g. including a variable time-step scheme). Moreover, the association and imple- 
mentation with existing standard, parallelization and multi-scale methods are also of 
interest. A demonstration code is available upon request to the author. 
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